function [MatDiag, SampleList, RKTree] = ConstructHMatrix(N, nLevel)

solvetime = 0;
RKmulttime = 0;

DIRECT = 1;
INDIRECT = 0;

nWidthTop = N / 4;
% nWidthBottom = 32;
% nLevel = round( log2( nWidthTop / nWidthBottom ) + 1);

maxRank = zeros(nLevel, 1);
maxRank(1) = 20;
maxRank(2:end) = 10;

nSample = maxRank + 10;
% SVDCut = 1d-6;
SVDCut = zeros(nLevel, 1);
SVDCut(1) = 1e-6;
for i = 2 : nLevel
  SVDCut(i) = SVDCut(i-1) * 2;
end

SampleList = ConstructSampleList( N, nLevel );
RKTree = ConstructRKTree(N, nLevel, SampleList);

display('Constructing RKTree...');
if(1)
  maxS = zeros(nLevel, 1);
  for iLevel = 1 : nLevel
    display( iLevel )
    tic
    
    
    nWidth = RKTree{iLevel}.nWidth;
    nGroup = SampleList{iLevel}.nSource / ...
      SampleList{iLevel}.nParallel;
    TempR = randn( N, N, nSample(iLevel) );
    TempQ = cell(RKTree{iLevel}.nRKBlock, 1);

    % Preparing random matrices
    for i = 1 : RKTree{iLevel}.nRKBlock
      TempQ{i} = zeros( nWidth, nWidth, nSample(iLevel) );
    end

    % Applying the kernel to the random matrices
    for iGroup = 1 : nGroup

      TempMatApply = zeros( N, N, nSample(iLevel) );
      for iParallel = 1 : SampleList{iLevel}.nParallel
	iSource = (iGroup-1) * SampleList{iLevel}.nParallel + ...
	  iParallel;
	Source = SampleList{iLevel}.Source( :, iSource );
	TempMatApply( 1+(Source(1)-1)*nWidth : ...
	    Source(1)*nWidth, ...
	    1+(Source(2)-1)*nWidth : ...
	    Source(2)*nWidth, : ) = ...
	    TempR( 1+(Source(1)-1)*nWidth : ...
	    Source(1)*nWidth, ...
	    1+(Source(2)-1)*nWidth : ...
	    Source(2)*nWidth, : );
      end
      TempMatRes = KernelApply( N, reshape(...
	TempMatApply, N * N, nSample(iLevel) ) );
      TempMatRes = reshape( TempMatRes, N, N, nSample(iLevel) ); 
      TempMatRKRes = RKTreeApply( RKTree, TempMatApply, N, iLevel, nSample(iLevel) );
      TempMatRes = TempMatRes - TempMatRKRes;
        
      % Extracting the blocks for each target and perform SVD
      for iParallel = 1 : SampleList{iLevel}.nParallel
	iSource = (iGroup-1) * SampleList{iLevel}.nParallel + ...
	  iParallel;
	for iTarget = 1 : SampleList{iLevel}.nTarget
	  iRKBlock = SampleList{iLevel}.RKIndex(iTarget, iSource );
	  Target = SampleList{iLevel}.SourceTarget( :, ...
	    iTarget, iSource );
	  TempT = TempMatRes(...
	    1+(Target(1)-1)*nWidth : ...
	    Target(1)*nWidth,  ...
	    1+(Target(2)-1)*nWidth : ...
	    Target(2)*nWidth, : );
	  [TempU, TempS, TempW ] = svd( ...
	    reshape( TempT, nWidth * nWidth, nSample(iLevel) ), 'econ' );
	  nRankU = min( max(find( abs(diag(TempS)) ./ max(abs(diag(TempS))) > SVDCut(iLevel) )), ...
	    maxRank(iLevel) );
	  if( nRankU == maxRank(iLevel)  )
	    disp('max rank reached')
	  end
	  if( SampleList{iLevel}.SymFlag(iTarget, iSource ) == DIRECT )
	    RKTree{iLevel}.RKBlockU{ iRKBlock } = TempU(:, 1:nRankU );
	    TempQ{iRKBlock} = TempT;
	  else
	    RKTree{iLevel}.RKBlockV{ iRKBlock } = TempU(:, 1:nRankU );
	  end
	  clear TempT TempU TempS TempW;
	end % iTarget
      end % iParallel	

      clear TempMatApply TempMatRes TempMatRKRes;

    end % iGroup

    % Calculating Sigma and form low rank approximation
    for iRKBlock = 1 : RKTree{iLevel}.nRKBlock
      Source = RKTree{iLevel}.RKSource(:, :, iRKBlock);
      Target = RKTree{iLevel}.RKTarget(:, :, iRKBlock);
      TempR1 = reshape( ...
	TempR( Source(1,1):Source(1,2), Source(2,1):Source(2,2), : ), ...
	nWidth * nWidth, nSample(iLevel) );
      TempR2 = reshape( ...
	TempR( Target(1,1):Target(1,2), Target(2,1):Target(2,2), : ), ...
	nWidth * nWidth, nSample(iLevel) );

      TempSigma = pinv( transpose(TempR2) * ...
	RKTree{iLevel}.RKBlockU{iRKBlock} ) * ...
	transpose(TempR2) * reshape( ...
	TempQ{iRKBlock}, nWidth * nWidth, nSample(iLevel) ) * ...
	pinv( transpose(RKTree{iLevel}.RKBlockV{iRKBlock}) * ...
	TempR1 );
      SSigma=svd(TempSigma);
      maxS(iLevel) = max(maxS(iLevel), max(SSigma));
      sizeSigma = size(TempSigma);
      if( sizeSigma(1) > sizeSigma(2) )
	RKTree{iLevel}.RKBlockU{iRKBlock} = ...
	  RKTree{iLevel}.RKBlockU{iRKBlock} * ...
	  TempSigma;
	RKTree{iLevel}.RKBlockV{iRKBlock} = ...
	  RKTree{iLevel}.RKBlockV{iRKBlock};
      else
	RKTree{iLevel}.RKBlockU{iRKBlock} = ...
	  RKTree{iLevel}.RKBlockU{iRKBlock};
	RKTree{iLevel}.RKBlockV{iRKBlock} = ...
	  RKTree{iLevel}.RKBlockV{iRKBlock} * ...
	  transpose(TempSigma);
      end
      clear TempR1 TempR2 TempSigma;
    end % iRKBlock
    
    clear TempR TempQ;

    maxS
    toc
  end % iLevel
end

% Extract diag part

display('Extracting diagonal...');
tic

% nWidth is doubled because we cannot sample two adjacent blocks.
nWidth = nWidth * 2;

MatDiag = zeros( N, N );
nParallel = SampleList{nLevel}.nParallel;

% This version is for the case when we have large memory.
if(1) 
  nSampleDiag = nWidth * nWidth;

  TempMatApply = zeros( N, N, nSampleDiag );
  for iSourceY = 1 : nWidth
    for iSourceX = 1 : nWidth
      SourceYInd = iSourceY : nWidth : N;
      SourceXInd = iSourceX : nWidth : N;
      ind = iSourceY + (iSourceX-1) * nWidth;
      TempMatApply(SourceYInd, SourceXInd, ind) = 1;
    end
  end

  TempMatRes = KernelApply( N, reshape(...
    TempMatApply, N * N, nSampleDiag) );
  TempMatRes = reshape( TempMatRes, N, N, nSampleDiag ); 
  TempMatRKRes = RKTreeApply( RKTree, TempMatApply, N, nLevel+1, nSampleDiag );
  TempMatRes = TempMatRes - TempMatRKRes;

  for iSourceY = 1 : nWidth
    for iSourceX = 1 : nWidth
      SourceYInd = iSourceY : nWidth : N;
      SourceXInd = iSourceX : nWidth : N;
      ind = iSourceY + (iSourceX-1) * nWidth;
      MatDiag(SourceYInd, SourceXInd) = ...
	TempMatRes(SourceYInd, SourceXInd, ind);
    end
  end

  clear TempMatApply TempMatRKRes TempMatRes;
end

% This version is when there is no such large memory.
if(0)
  nSampleDiag = nWidth;

  for iSourceX = 1 : nWidth
    TempMatApply = zeros( N, N, nSampleDiag );
    for iSourceY = 1 : nWidth
      SourceYInd = iSourceY : nWidth : N;
      SourceXInd = iSourceX : nWidth : N;
      ind = iSourceY;
      TempMatApply(SourceYInd, SourceXInd, ind) = 1;
    end
    TempMatRes = KernelApply( N, reshape(...
      TempMatApply, N * N, nSampleDiag) );
    TempMatRes = reshape( TempMatRes, N, N, nSampleDiag ); 
    TempMatRKRes = RKTreeApply( RKTree, TempMatApply, N, nLevel+1, nSampleDiag );
    TempMatRes = TempMatRes - TempMatRKRes;
  
    for iSourceY = 1 : nWidth
      SourceYInd = iSourceY : nWidth : N;
      SourceXInd = iSourceX : nWidth : N;
      ind = iSourceY;
      MatDiag(SourceYInd, SourceXInd) = ...
	TempMatRes(SourceYInd, SourceXInd, ind);
    end
    clear TempMatApply TempMatRKRes TempMatRes;

  end

end

toc

